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Abstract . In all three principal methods of measuring the earth’s gravitational 
field — gravimetry, satellite orbit perturbations, and astrogeodetic networks — both mathe- 
matical theory and observational techniques have been developed in recent years to more 
than sufficient accuracy to define and determine the irregularities of the field significant 
on a planetary scale. The common defect of all three methods is inadequate distribution 
of observations. At present, the different methods agree only for major features such as 
the Indian Ocean minimum, and the normalized spherical harmonic coefficients of the 
potential C nm , S nm are known to be of an order of magnitude of about ±1.3 X 10 -6 /(7i — 1) 
with an uncertainty of about ±0.6 X 10“ 8 for the lesser al terms. Gradual improvement 
is anticipated with additional observations and the application of better statistical tech- 
niques made possible by modern computers. 




INTRODUCTION 

The subject of this review is the determination of the earth’s gravitational 
field, given measurements that are adequate in accuracy but incomplete in their 
distribution. The emphasis will be on the wide-scale variations of the field rather 
than on local variations. In accordance with this emphasis, we shall usually employ 
the spectral representation of the gravitational field rather than the spatial. 

Detailed discussion will be attempted only for developments since the advent 
of artificial satellites in 1957. Earlier developments are described adequately in 
publications generally available: the texts by Heiskanen and Vening-Meinesz [1958], 
Jeffreys [1959], and Bomford [1962]; the Handbuch der Physik articles by Garland 
[1956] and Jung [1956]; and, recently translated into English, the important 
treatise by Molodenskiy , Yeremeyev , and Yurkina [I960]. 



MATHEMATICAL EXPRESSION OF THE GRAVITY FIELD 

The reference figure. The physically logical reference figure is that of a rotating 
fluid in equilibrium. The principal variable involved in such a model is the radial 
density distribution. The theoretical investigation of rotating fluid models of the 
earth is still being pursued, principally by Leder sieger [1962; see this reference for 
earlier papers]. 

The more mathematically tractable ellipsoid of revolution, which differs by 
quantities of the order of 10~ 6 from a rotating fluid, is generally used in practice 
for a reference figure. A closed expression for the gravity (i.e., gravitational attrac- 
tion plus centrifugal force), first obtained by Pizzetti [1894], requires the use of 
ellipsoidal harmonics. For practical application, the present accuracy of determina- 
tion of the gravity field makes it more convenient to use series developments in 
spherical harmonics. For an ellipsoidal figure specified by an equatorial semimajor 
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axis a, an equatorial gravity y e> a polar semiaxis 5, and a rate of rotation w, formulas 
of sufficient precision for the subsequent discussions in this review are [Jung, 1956, 
pp. 542-562; Heiskanen and Vening-Meinesz, 1958, pp. 51-53; Bomjord, 1962, 
pp. 414-416] 



/ = (a - 

6)/a 


(1) 


m = co 2 a/y. 


(2) 

kM 

= a\,[l -/+§«-** mf + 0(f)] 

(3 ) 

(1 m — 

/ - tfmf) sin 2 *„ 

+ (if - 

- fmf) sin 2 2*, + 0(f)] 

(4) 

e = afl 

- (/ + if) sin 2 *, 

, + §f sin 4 *, + 0(f)] 

(5) 

J 2 = 1 

1/(1 -if)- Mi 

— § m - 

- 1/] + 0(f) 

(6) 


* 

ii 

I 

> 

' — K 

1 

5m) + 

0(f) 

(7) 

kM 

- P 2 ( sin *) 


) P 4 (sin *) + 0(f)] 

(8) 


where k is the gravitational constant, M is the mass, <j> is the geocentric latitude, 
4> 9 is the geodetic latitude, r e is the radial spherical coordinate of the surface of 
the ellipsoid, U is the gravitational potential, and the P n are Legendre poly- 
nomials. Equations 3 through 7 are essentially those given by Helmert [1884] and 
are used by most geodesists. Another set derived by DeSitter [1938] is preferred 
by most astronomers. 

The expression of the external field of an ellipsoid of revolution has recently 
been developed in series to higher-order terms by Cook [1959], Hirvonen [I960], 
and Lambert [1961] and expressed in closed formulas by Caputo [1963]. 

Variations in the field : first approximation. To accommodate departures of 
the gravitational potential of the earth from the reference figure, the field can be 
expressed as a sum of spherical harmonics: 

K= — [l+ 2]f-) 2 -P»m(sin <t>){C nm cos mX + §„„ sin m\} 1 (9) 

' L «- 2 W m- 0 J 


The gravitational potential V is positive, following the sign convention of 
astronomy and geodesy. The absence of n = 1 terms from (9) means that the co- 
ordinate origin coincides with the center of mass. 

P nm (sin <j>) is the Legendre associated polynomial. The overbar signifies a 
normalization; the one we shall apply in this review makes 

J J j~P nTn (sin J cos <M<M\ = 4?r (10) 
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where 5 0m is the Kronecker delta, and k is the integer part of (n — m)/ 2. 

The coefficients C n „, S„ m are independent parameters 

C™ = -J 2 /VI « -0.0004842 (12) 

All other G nm , S n „ are 0(10 6 ). The magnitude of / (or J 2 or C 20 ) means that varia- 
tions of the field referred to the ellipsoidal surface r e expressed by (5) can be ex- 
pressed by the C nm , S nm with a proportionate error of only H)~ 3 . The variation most 
commonly referred to the ellipsoid is the geoid, i.e., the equipotential surface of 
the family defined by (9) for which V is equal to the U of (8) on the ellipsoidal 
surface in (5). Expressed in this manner, the geoid height is 

N = (V - U)/y = T/y (13) 

The ellipsoidal surface r e is an equipotential for the gravitational potential U 
plus the centrifugal force pote ntial \p : 

. 2 

~ ~ r 2 cos 2 4> (14) 

Measurements of the gravity acceleration can be referred to the geoid through 
leveling observations. To compare the observed acceleration g with the accelera- 
tion y on the reference model, we must therefore allow for the height difference N: 

Ag = g(r* + N) — y(r e ) 

= -[^1 .+[**&**+**£* L < i5 > 

where the derivatives are normal to the ellipsoid. 

Using (13) and taking the spherical approximation —2y /r for dy/dr } we obtain 

Ag = (dT/dn) — 2T/r (16) 

** * 

The last term in.(lfi) is usually referred to as Bruns’ term. 

. t For the coefficient A nm , B nm of a particular term in the spherical harmonic 
expression of A g, we obtain from (9), (15), and (16) 

A nm ,B nm = y(n - l){C nm , S nm } (17) 

■4.-.,.- To relate the geoid height A 7 at a particular point to the anomalies A g over the 
surface of the earth, we develop the anomalies in spherical harmonics about the 
poi^t as a pole and use the orthogonality relationships between spherical harmonics 
to obtain: • 

r r 

‘ N = ~ / S( cos 6) A g da (18) 

TrTTy J sphere 

where 

S( cos e) = E p,( cos e) (19) 

n-2 n — X 

P n (cos 6) is the normalized zonal harmonic. The closed form of S (cos 0) is Stokes’ 

function [ Stokes , 1849]: 

6 ( 6 8 
S ( cos 8) = esc - — 3 cos 8 In \sm - + sin 2 -J — 6 sin - + 1 — 5 cos 8 (20) 
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The similar function for the slope of the geoid, or deflection of the vertical, was ob- 
tained by Vening-Meinesz [1928] and de Graaff-Hunter [1935]. 

If the coordinate system were rotated and the generalized additional theorem 
applied to the spherical harmonics, the coefficients C nm , S nm referred to the new 
axes would be expressed entirely as functions of the coefficients of degree n referred 
to the old axes. Furthermore, the quantity 

= £ { C n J + S n J} (21) 

m 

remains invariant under rotation [Kaula, 19595]. Thus, the ‘degree variance’ <r„ 2 
is a measure of the amount of variability in a certain wavelength that is inde- 
pendent of the coordinate axes’ directions. By summary a 2 over n from 2 to in- 
finity, we obtain the mean square amplitude. 

The foregoing relationships were all derived on the assumption that the func- 
tions were harmonic, which requires that the geoid be external to all masses. Since 
the actual geoid is not external to all masses, hypothetical transfers of mass must 
be made before ‘reduction to sea level 7 of gravity observations and application of 
Stokes’ function. This problem has historically occupied most of the attention of 
physical geodesy [Heiskanen and Vening-Meinesz , 1958, pp. 147-186; Bomfordj 
1962, pp. 416-446]. Advocates can still be found for most of the reductions that 
have been proposed in the past. The free air anomaly, which implies that the ex- 
ternal topography is condensed in a thin layer just below the geoid; smoothed 
modifications of the free air anomaly [de Graaff-Hunter , I960]; the isostatic anomaly, 
which implies that topographic excesses and deficiencies are redistributed uniformly 
through the crust (Pratt-Hayford compensation) or transferred to the crust-mantle 
interface (Airy-Heiskanen compensation); and even the Rudzki anomaly, which 
implies that the external masses are replaced by internal masses yielding the same 
geoid [ Tengstrom , 1962] are all still applied for geodetic purposes, while the Bouguer 
anomaly, which implies complete removal of the mass excess and deficiencies, is 
most often applied in geologic interpretation. 

The free air anomaly is the most often used, mainly for the practical reason 
that it does not require detailed evaluation of the topography (except for the local 
area). As Levallois [1962] emphasized, all the methods of reduction can be con- 
sidered as equivalent provided that the ‘indirect effect’ (the effect on A g of the 
geoid shift that occurs because of the implied mass transfer) of each method is 
properly taken into account. Hence the practical distinction becomes which type 
of anomaly can be considered as most ‘representative’; i.e., for which type can the 
observed point values be interpolated most easily. 

Isostatic anomalies are generally regarded as the most representative, but 
they are also the most laborious, since it is difficult to automate estimation of topo- 
graphic elevations. Work is continuing to facilitate Airy-Heiskanen isostatic re- 
ductions, such as, e.g., the reduction maps for the distant zones of the Isostatic 
Institute [Karki et al., 1961]. For analytic application of the Airy-Heiskanen iso- 
static reduction to spherical harmonics, formulas of Jung [1952] are convenient. 
The isostatic reduction can be expressed as a power series in h(<p, X), the height of 
the solid surface, and X), the depth of the water. By expressing w as a surface 
coating, Jung develops the isostatic correction in terms of the spherical harmonic 
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expansion of the topography plus water to quadratic terms in the topography 
of the form: 

oo 

A 9i. = Z \e»[El + UE \ } (22) 

»*»0 

where [E] n and \E 2 ] n are the sums of terms of degree n in the harmonic expression of 

E = h- (1.03/<r o)w (23) 

and its square, and e n and / n are functions of the depth of compensation and the 
crustal and mantle densities. The functions e n and f n contain terms of the form 
{1 — (1 — i) 2n+1 }, where t is the ratio of the compensation depth to the earth's 
radius, and so (22) is probably imprecise for higher values of n for this reason as 
well as because of the lack of higher powers of E. For the purposes of this review, 
(22) can be considered a means of estimating &g where E is known but not g. 

Variations in the field : higher approximations . As in any problem, further 
elaboration can be made either in the reference model or in the expressions of 
departures from the model. The first type of elaboration was made by Sagrebin 
[1956] who refined Stokes' formula (equation 18) to an accuracy of order f by using 
elipsoidal harmonics; his solution has since been improved and corrected by 
Molodenskiy et aL [1960] and by Bjerhammar [1962b]. These investigators found 
that the error in Stokes' formula was indeed of the order of /, i.e., 20 cm or less in N. 

The principal theoretical developments in physical geodesy in recent years 
have been in nonlinear expression of the departures from the model by considering 
the problem at the physical surface of the earth rather than at the geoid, thus 
eliminating the inaccuracies due to ignorance of the density within the earth. 
These developments were initiated in a paper by Jeffreys [1931], and have been 
carried forward principally by Molodenskiy [1945, 1948]. Full details are given by 
Molodenskiy et ah [1960], and a brief summary by Molodenskiy [1962]. The deriva- 
tion of the method currently applied in the USSR starts with the expression of the 
disturbed potential T in the form of a surface layer <p on the physical surface of the 
earth. Taking the derivative of T normal to the physical surface, using the Bruns' 
equation 16, and allowing for the slope of the surface with respect to the equi- 
potential, we obtain an expression for the free air anomaly Ag F in terms of integrals 
over the equipotential surface S. When we assume that the projection of the surface 
layer <p onto a sphere S, of radius R, can be expressed in terms of a power series 
over the discrepancy between S and S , the solution is obtained as a successive 
approximation in four cycles, n = 0 through 3, over (setting negatively subscripted 
variables zero) 

Gn = R J Z=_^> ^ Xn _ 1 _ 2 (h _ h 0 )Xn-2 ~ d<J 

+ 2 ttx „~2 tan 2 a + 8 0n Agr (24) 

x - = S + (ifr / G »w cos *> - « d(J ( 25 ) 

where h 0 is the ‘normal' height at the point where the quantity on the left is com- 
puted, h is the ‘normal' height at the location of the quantities within the integral, 
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r 0 is the chord distance of the sphere between the two points, a is the slope of the 
surface at the point of computation, and the integration is over the unit sphere. 
The ‘normal' height differs from the orthometric height in that the normal gravity 
is used in place of the true gravity g. Then to obtain T , the disturbance of the 
potential: 

Tn== &f G "^ cos 0) - il - y / - - 7/°" x»-s da ( 26 ) 

T = £ T n (27) 

n 

Similar formulas are derived for the deflection of the vertical. 

Theoretical developments similar to Molodenskiy's have been made by 
Levallois [1957], Arnold [1959a, 5], de Graaff-Hunter [1960], Hirvonen [1960], Moritz 
[1961], Bjerhammar [1962a], and Cook [19636]. Summaries and comparisons are 
given by Arnold [1960a], Hirvonen [1961], Tengstrom [1961], and Molodenskiy 
et aL [1962], who criticize the approximations made in some of the other theories. 

Examples of the numerical significance of the improved theory are given by 
Molodenskiy et al [1960], Arnold [19606], Tengstrom [1961], and Pellinen [1962]. 
The conclusion is that the difference in geoid height is always small, of the order 
of tens of centimeters at most, but that the effect on the deflection of the vertical 
can be appreciable, several seconds of arc, in extreme situations. Of particular 
interest for the determination of the long-wave components of the field is the com- 
putation by Pellinen [1962] of the effect of the G x term from (24) on the low-degree 
harmonics. He found that the rms effect up through degree 3 was d=0.20 mgal 
on the normalized A g coefficients, or about d=0.15 X 10 -6 for the potential coeffi- 
cients C nmj S nm . 


STATISTICAL CONSIDERATIONS 

Random processes referred to two types of continuums are of concern: time 
series, which comprise the perturbations of a satellite orbit, and distributions over 
a spherical surface, which comprise the variations of the earth's gravity field. 

We shall first discuss heuristically the general questions of spectral analysis 
and estimation by quadratic sum minimization, and then apply the conclusions 
to the particular cases of time series and distributions over a spherical surface. 

Spectral analysis and quadratic sum minimization . Discrepancies f(s) of ob- 
servations from a mathematical model that are small enough to be considered as 
linear can be represented as a linear transformation of a set of parameters x: 

f(s) = c T (s)x s e S (28) 

The superscript T denotes the transpose. The coordinates s may be of any 
number of dimensions. S is a subspace, not necessarily connected, of a total space T. 
Some of, or all, the parameters x are members of an enumerable infinite set which 
are orthogonal over T, i.e., the square array 
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is diagonal (dr is an element of volume in T). The functions c(0 are usually nor- 
malized so that the diagonal elements of A are either 1 or a constant times J T dr . 
x can be divided into subsets x n such that c n T (t)x n is invariant under rotation of the 
coordinate systems. The quadratic sum 

* n 2 = X n T Xn (30) 

over one of these subsets is the degree variance, already mentioned in (21) for the 
case of T which is a spherical surface. For many phenomena, only the <r„ 2 are of 
interest, and not the individual components of x n . Usually it is assumed that the 
random process is isotropic, i.e. the covariance 

K(r,s) =E{f(r)Ks)\ (31) 

where E denotes the mean value over T and is a function of only the distance (the 
length along the geodesic) between r and s. Under isotropy, the covariance can be 
expressed as 

K(r, s) = Z KcMr.W (32) 

n 


where d r . is the distance between r and s f c n0 is a member of c n , and k n is constant. 
The process of spectral analysis is that of forming numerical estimates of if(r, s) 
and then using the orthogonality expressed by (29) to determine the a n 2 . The 
principal problem in this process is sampling: (1) the distribution of sample points, 
considering the wavelengths corresponding to the highest degree a 2 anticipated 
to be significant (the ‘aliasing’ problem); and (2) the weighting, if any, applied 
to take into account the incomplete extent of the observed subspace S with respect 
to the total space T (the ‘window’ problem). 

If the individual components of x n are wanted, or if some of the parameters 
of x are not coefficients of orthogonal functions, then spectral analysis is insufficient 
and quadratic sum minimization must be applied. In this case, the estimates 
K(r , s) are still made, but the x are determined so that, subject to (28), 



x r c (r)K \r, s)c T (s)x dr ds — minimum 


where the inverse K *(r, s) is a function such that 


(33) 


f K7\r, s)K(s , t) <k = $(r, 0 
Js 

where 5(r, t) is the Dirac delta function. The solution is 


(34) 


x = W f f c(r)ir x (r, *)/(«) drds (35) 

Js Js 

where W is a symmetric matrix whose elements are the estimated variances and 
covariances of the x’s: e.g., <r n 2 /AT for the orthogonal function coefficients, where 
N n is the number of terms of degree n. 

The minimization of (33) is most generally proved as the determination of the 
projection c T x in the subspace specified by the chosen elements of x of the vector 
/(r) in a Hilbert space with metric fU^r, s), or measure K(r, s). 
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More familiar is the form of (28) and (33) where (1) /(«) is evaluated at a 
finite set of points, and (2) for some of the parameters in x there is no a priori 
estimate of the variances. Also the covariance matrix W of the x may be known 
rather than K, the matrix equivalent of K(r, s). In this case, the set of points can 
be expressed as a vector, and the varianceless part of x broken off as a separate 


vector z: 

Cx + Mz = f (36) 

(f — Mz) r K _1 (f — Mz) = x T W _1 x = minimum (37) 

At least one element is nonzero in each row of C. 

Solution by the method of Lagrangian multipliers obtains: 

x = WC T K \I - M(M r K“ 1 M)‘ 1 M T K“ 1 )f (38) 

z = (M r K“ 1 M)" 1 M r ET 1 f (39) 

where 

K - CWC r (40) 

and I is the identity matrix. The covariance matrix of z is 

V = (M r ET l M r 1 (41) 

and of the corrected x 

U = W - WC r K _1 (I - M(M T K _l M)“ 1 M r K" 1 )CW (42) 


The rows of WC r (and columns of CW) in (38) and (42) need be only those of the 
elements of x of interest. If K is calculated by (40), however, the W must include 
all parameters of appreciable effect. Conventionally, the elements of x are some- 
times referred to as ‘observations' and the elements of z as ‘parameters' ; mathe- 
matically, the distinction is whether or not covariance W can be preassigned. 
Sometimes there are added to the system of (36) and (37) ‘side conditions': 

Nz = k (43) 

is a minor complication that can be removed by eliminating a corresponding num- 
ber of elements in z. 

Ordinary least squares is the case in which 

C - W = I (44) 

Prediction and interpolation by linear regression is the case in which it is 
desired to estimate / in the subspace T — S: 

/(r) = c r (r)x rzT - S (45) 

or, letting g be a set of values of / in T — S and D be the corresponding array of 
c r 's, and using (38) (assuming no z) 

E\ g} = DWC r K‘ 1 f = Bf (46) 

A result that is also arrived at through solving the Wiener-Hopf equations for 

the regression coefficients B is 


DWC r = BCWC r 
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or 


K = BK„ (47) 

Sometimes the solution for x or z is made in stages, with different condition 
equations 36 in each stage, either because the nature of the process is evolutionary 
(e.g., orbit predictions) or because a single stage computation would be too un- 
wieldy. In this case, the vector x for a particular stage can be considered as con- 
sisting of three parts: (1) new observations x XJ with associated coefficient matrix 
C x and covariance matrix W x ; (2) old observations x„, with associated coefficient 
matrix C„ and covariance matrix W„: 


W y = P v UoP/ (48) 

where U 0 is the U computed by (42) in the previous stage, or a submatrix thereof, 
and P„ is a propagation or transition matrix taking into account change of the 
parameters (or ‘state variables’) specifying the reference model; and (3) previously 
estimated parameters x„ with associated coefficient matrix C, and covariance 
matrix W, : 

W t = P.VoP/ (49) 

where V 0 is the V of (41) or a submatrix thereof. The C„ C v , C, are combined to 
form a new C and the W„ W v , W, are combined to form a new W for another 
solution according to (36) through (41). In calculating the residuals f for each 
new stage, the corrections x, z from previous stages are incorporated in the reference 
model. 

The foregoing discussion is a consequence of attempts to combine ideas of 
time series analysis of Bartlett [1956] and Parzen [1961] with least squares and its 
generalization as given by Arley and Buck [1950] and Brown [1955, 1957]. The 
results can probably be found to be subsumed in the general treatments of Bochner 
[1955] and Yaglom [1961] by those able to cope with the recondite mathematics 
therein. The derivation of (36) through (42) has recently been discussed in detail 
by Steam and Richardson [1962]. The linear regression prediction of (45) through 
(47) is applied by Moritz [19625] to observations on a plane and by Parzen [1961] 
to time series in continuous, rather than discrete, form. Special cases of the staged 
computation of x x , x y , x z with W„, W f calculated as in (48) and (49) are: (1) the 
optimal prediction of Kalman [I960]: W„ C x , and M are 0; and (2) the Bayes 
estimation of Parzen [1962a] and the preassigned covariance of Kaula [1961c]: 
W v and C y are 0; P, is I; and C x , C„ M, f have the forms 

C.-$ C, = {i} M = {f} 

Another variation of estimation by quadratic sum minimization is Gram- 
Schmidt orthogonalization [e.g., Robinson , 1959]. The function f(s) of (28) is 
represented as 



f(s) = d r y stS 


(51) 
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where the functions d are orthogonal over S, not T; i.e., 


d(s)cT(s) dr = I 
Js 

(52) 

the identity matrix. Setting 


Bd = c 

(53) 

we get 


BB T = [ cc T dr 
Js 

(54) 


B is thus not unique. The Gram-Schmidt orthogonalization process finds a unique 
B by requiring that, for the arbitrarily selected sequence of functions in c, B is 
triangular. Letting the subscripts on d i} c { denote ordering in the arbitrary sequence 
Gram-Schmidt orthogonalization yields 

di(s) = c x (s) (55) 

*-i d k (s) 1 d k c k dr 

di(s) = Ci(s) - £ ^ i > 1 (56) 

*"* f d k 2 dr 

If the x of (28) has a finite number of members, x can be determined by simple 
least squares; i.e., in (39) z represents x, K is I, and M is the array of values of c 
and f is the array of values of / for a set of points or volume element means. In 
practice, x has an infinite number of members, and it is not known beforehand 
which are significant; also, the disposition of S within T may be such that M r M 
approaches singularity, so that it is computationally impracticable to invert. If 
Gram-Schmidt orthogonalization is applied, then by (52) M r M is diagonal and 
(39) degenerates to a simple Fourier analysis in which the coefficients are deter- 
mined one by one until it is decided that nothing more of significance is being found. 
Gram-Schmidt orthogonalization thus is a systematic method of estimation that 
avoids the inversion (M^M)" 1 of simple least squares or the inversion K" 1 of 
quadratic sum minimization (equation 35 or 38 with M = 0). It does not, however, 
abolish the ‘window’ effect — the distortion of estimates of lower-degree coefficients 
by higher-degree variations — and by avoiding the assumptions of stationarity and 
isotropy of spectral analysis and quadratic-sum minimization it is accordingly a 
less effective predictor of the variations in the subspace T — S. 

A technique similar to Gram-Schmidt orthogonalization in principle is step- 
wise least squares : the parameters in x are determined one each or one subset at a 
time, and the contribution of a particular set is subtracted from f before determining 
those of the next set. This method is commonly employed in excessively large 
problems when the parameters are not coefficients of orthogonal functions. It 
could perhaps be considered as a special case of the staged least squares (equation 
50) in which W z is set zero, but differs in that the same condition equations 36 
are used repeatedly with different sections of the M matrix set zero, which implies 
that the situation must approach orthogonality for the procedure to be justifiable. 
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Time series. Probability models of time series, and the analysis of continuous 
or uniformly spaced observations thereof, have been developed in considerable 
detail [Bartlett, 1956; Grenander and Rosenblatt , 1957; Blackman and Tukey, 1959; 
Parzen , 1961]. Applications to geophysics have been discussed by Holloway [1958], 
Munk and MacDonald [1960], and Van Isacker [1961]. 

The time series of significance in the present problem are vector quantities 
characterized by a discrete spectrum of fixed phase and amplitude plus a continuous 
spectrum of noise [Kaula, 1963a]: 


x f (t) = a fn exp [i\ n t] + J 0 f (oS) exp [uoZ] da>j 


(57) 


where the subscript / denotes a vector comp onen t and n denotes one of a set of 
frequencies; (R denotes the real part, i is \/ — 1, and the vectors a fn , £,(c c) are 
complex. The X„ are all known; the a fn are linear functions of a set of parameters p h 
less than N in number, which it is our problem to determine. The duration T of 
observation of the time series is such that the discrete spectrum stands out above 
the noise; i.e., 


/3 / (co)(2x/oj7 7 ) « a fn co > X» (58) 

but 

j3 / (aj)(2x/ajT') a fn some a> <<C X n (59) 

If the x f (t) were continuously observed with observational errors random, the 
condition expressed by (58) would assure that the p t could be accurately deter- 
mined. In the problem of interest, however, the x f (t) are incompletely and inter- 
mittently observed, and furthermore the observations are affected by unknown 
biases; i.e., x f (t) must be replaced by 

*,m - no g; [„,(» - ,.] (6o> 

where I(t) is unity during observation and zero at other times, and the number of 
components denoted by j is less than those denoted by /. Under these circumstances, 
the effective orthogonality expressed by (58) is destroyed and the condition ex- 
pressed by (59) greatly increases the length of record required to estimate the 
parameters p h q m by ordinary methods such as simple least squares, as well as 
increasing the chance that sources of systematic distortion could be hidden in the 
noise >). 

The problem of spectral analysis of time series with missing observations has 
been treated by Jones [1962] and Parzen [19626]. These discussions treat I(t) as an 
‘amplification factor': 

x(t) - I{t)y{t) (61) 

and show that for covariance estimates for any lag v 

ft(p) - RM/Ri(v) (62) 

The condition expressed by (59) suggests, however, that the time series itself 
would be useful to determine only the covariance due to the low-frequency con- 
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tiuuous part of the spectrum, and that since the da fn /dpi are known it would be 
better to obtain estimates of the high-frequency discrete spectrum effects from 
independent estimates of the variances of the p ( ’s. The practical problem then be- 
comes the large dimension of the covariance matrix which must be manipulated 
in computation. 

Distributions over a spherical surface. Spectral representation on a spherical 
surface is discussed by Schoenberg [1942], Obukhov [1947], Kaula [19596], and J ones 
[1963]. In this case, the coordinates s of (28) are latitude and longitude (<f>, X) and 
the functions c T (s) are surface spherical harmonics. With the normalization specified 
by (10) and (11), and the assumption of isotropy, the covariance of (32) becomes 

K(r,s) = £ P„,(cos 6„) (63) 

n v2 n + 1 

The distributions over a sphere significant in the present problem are scalar 
quantities characterized by spectrums with significant contributions from wave 
numbers as high as 4000. Hence the sample to determine the covariances K(0 ra ), 
and thence to estimate the cr rt 2, s for small n* s should be large and randomly distri- 
buted to avoid aliasing [Shapiro and Silverman , I960]. Furthermore, in (38) and 
(39) the matrix K employed must be calculated directly from the estimated co- 
variances K(0 r ,) and not from W, as stated in (40). 

The cross covariance between two different variables /, h on a spherical surface, 
assuming isotropy, will be 

K,M = E -W==T ^( cos O (64) 

n v2 n + 1 


where 

MW I < (65) 


If / is known over part (S) of the sphere, and h is known over all ( T ) the sphere, 
the best estimate of a coefficient C nm (f) will be the weighted mean of the estimate 
obtained from the / in S by (35) and the estimate obtained using C nm (h) and 
the cross-degree variance <r n (jh ) : 


E[C nm (f)} 

1 r <r n 2 (f) 
Ps + Ph L Ps 2n + 1 


u 


c nm {r)K-\r, 8)1(8) dr ds + Pk ^ <?--W 


\h) 


(66) 


P* rt 


<r 2 {C nm (J) \ ftS] 


Jt+i {' - I, L s)c - is) * * 


(67) 


1 1 _ 

PK * 2 {CUf) I C n Jh)\ 1 / _ M6)] 2 \l 

L^n r " (f) _ a n \h )} J 


( 68 ) 
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Similarly, for the optimum estimate of / at a set of points in T — S a modifica- 
tion of (46), (47) can be used: 

E\g) = [V ,,- 1 + V.rT 1 [V f/ - , X,,K,r 1 f + V/K.Ai-’h] (69) 

where 

V., = U„ - (70) 

= U„ K^Eu (71) 

in which 

u„ = I* 2 {/} (72) 

If it is assumed that / can be estimated from h by some physical rule such 
as the isostatic rule of (22), then the foregoing statistical treatment should be 
applied to / less that part of / accounted for by the physical rule. This treatment 
is to be emphasized if the parameters of the rule were determined from samples of 
extreme, rather than average, characteristics, as were those of the isostatic reduc- 
tions [Heiskanen and Vening-Meinesz, 1958, pp. 187-2211. 

Spherical harmonic analysis is most common in geomagnetism, where the 
standard methods have been to interpolate to a latitude and longitude grid, to 
fit Fourier series to parallels, and then to fit Legendre functions to the parallel 
coefficients [< Chapman and Bartels , 1940]. Fougere [1963] determines polynomials 
by Gram-Schmidt orthogonalization and then determines the spherical harmonics 
corresponding to those polynomials found significant. The methods applied in 
geomagnetism, bypassing the covariance matrix, seem possible, however, only 
because of the relatively small contribution of higher harmonics. 

USE OF GRAVIMETRY 

Observing system . A review of gravimetric techniques will appear shortly 
[LaCoste and Harrison , 1964]; see also Harrison [1962]. Of principal concern in the 
present review are the distribution of gravimetric data and the magnitude of 
possible systematic errors. 

The distribution of available gravimetry is shown in Figure 1, which is based 
on the map of Uotila [1961] updated by information from Worzel et al. [1963] 
and McCahan [1963]. The principal defect of the distribution has always been, 
of course, the correlation with topography, both on the wide scale shown by Figure 1 
and locally, as, for example, in observations being conducted mainly in the valleys of 
mountainous areas and on the islands in archipelagoes. The principal difficulty in 
improving worldwide coverage is the development of economic airborne gravimeters 
with accurate velocity determination or seaborne gravimeters capable of obtaining 
accurate results in a reasonable rough sea. 

Systematic errors possibly affecting series of observations which either are 
the principal data for an appreciable area or provide the reference or calibration 
net for other observations include: (1) connection of national and other gravi- 
metric systems to the world reference network; (2) calibration: the scale factors 
utilized for gravimeters; (3) referral of seaborne measurements to land reference 
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stations; (4) navigation errors affecting the position and Eotvos correction for air- 
borne and seaborne gravimetry; and (5) horizontal acceleration and level error 
effects on airborne and seaborne gravimetry. Errors in categories (1), (2), and (3) 
are discussed by Woollard [1961; Woollard and Rose, 1963]. Woollard [1961] concludes 
that the principal gravimetric systems are connected with an average error of 
about ±0.4 mgal on the basis of gravimeter connections since 1948 and pendul um 
connections since 1953. Morelli [1963] is more pessimistic, citing as an example a 
range of 3.2 mgal in values obtained for the Rome reference station by various 
observers since 1960. The principal improvement needed is more reliable north- 
south pendulum lines to calibrate gravimeters. 

The referral of measurements at sea to land reference systems appears to be 
poorer: Woollard [1961] concludes that Vening-Meinesz’s pendulum values are 
probably reliable to ±2 mgal through measurements made on land with the 
pendulums, even though 13 of 24 comparisons with Vening-Meinesz’s harbor 
measurements differ by 5 or more mgal. Comparisons from track crossings of dif- 
ferent cruises are somewhat better, usually averaging ±5 mgal, which includes 
other effects. Navigation error is generally considered to be the principal source 
of error in air and sea measurements, limiting their random accuracy to ±3 to 
±5 mgal [Harrison, 1962]. There does not seem to be any estimate of systematic 
navigation error. Horizontal acceleration and level errors dependent on the heading 
of the ship relative to the sea of several milligals have been reported [Allan et al, 
1962]. Serious systematic effects of this sort, however, must be less than the ±5 
mgal of crossing comparisons with submarine measurements avoiding sea swell 
effects. 

It appears that the most likely source of significant systematic error in gravi- 
metry is in the connection to the reference system of a cruise of marine measure- 
ments which constitute the sole data in an appreciable area, but that this error is 
less than 5 mgal. 

The absolute value of gravity is needed for comparison of results obtained 
from gravimetry with those from satellite orbits. The correction to the Potsdam 
standard is currently estimated to be about — 13 ± 1 mgal. Improvement is antici- 
pated from several determinations in progress [Cook, 1963a]. 

An important practical task is the collection and processing of the millions 
of gravimetric observations that have been made. The leader in accomplishing 
this task has been W. A. Heiskanen, who established the gravity centers at the 
Isostatic Institute in Helsinki and the Ohio State University in Columbus. De- 
scriptions of the collection and processing are given by Uotila [1960] and Heiskanen 
[1962]. Other collections of gravity data exist at the U. S. Naval Oceanographic 
Office in Washington and the Bureau Gravimetrique Internationale in Paris. 
Standardization of punched card format for single observations is being worked 
out by these centers. To exploit fully the existing gravimetry for determination 
of the worldwide field there is still needed agreement on a standardization of the 
local treatment to obtain area means for, say, 1° by 1° or 100-km by 100-km 
blocks, as described in the next section, and of the recording thereof on punched 
cards or other automated form. 

Local treatment . Gravity anomalies calculated at observation points cannot 
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be used directly in the determination of covariances, harmonic coefficients, and 
geoid heights as described in preceding sections of this review because of the needs 
to smooth out the high amount of local variability; to remove, as much as possible, 
the effect of correlation of observation distribution with local topography; and to 
keep the computations to a manageable size. The first two of these desiderata indi- 
cate the use of a gravity reduction which takes into account the topography. 
However, only Tanni [1948] and Uotila [1962] have utilized isostatic anomalies. 
All other calculations of the geoid from gravimetry have placed greater emphasis 
on the third of the desiderata and have applied at most a simple linear correlation 
formula to reduce free air gravity anomalies to the mean elevation of the area 
they represent : 

Ag m — a Q bh (73) 

where h is the mean elevation. Hirvonen [1934] and Dubovskiy [ Molodenskiy , 1945] 
did not apply any such correction in their early determinations; Jeffreys [1943] and 
Zhonglovich [1952] applied the correction to the anomalies for 10° by 10° squares, 
which were the direct mean of observed values in 1° by 1° squares; Heiskanen 
[1957], Kaula [1959a, 6], and Uotila [1962] used mean anomalies of 1° by 1° squares 
corrected to mean elevation and then applied step-by-step, or Markov, extrap- 
olation and interpolation to obtain the mean values of 5° by 5° squares. The 
formation of the mean anomalies for 1° by 1° squares was done largely by Uotila 
[1960], who found mean values of b of 0.118 mgal/meter on land and 0.069 mgal/ 
meter at sea. In cases where the distribution of observations with respect to 
elevation in a 1° by 1° square was insufficient to determine a gradient, Uotila used 
the gradients from nearby squares or the mean gradients of all squares. 

Calculations of isostatic anomalies for much of the available gravimetry have 
been made at the Ohio State University, the Isostatic Institute [ Heiskanen , 1962], 
and the Bureau Gravimetrique Internationale, but published maps of the reduced 
anomalies have been limited to a few areas such as Europe and North Africa 
[Coron, 1962]. 

The development of computers encourages the application of analytic tech- 
niques to replace graphical methods for smoothing, averaging, etc. There have 
been several papers on use of Fourier series [Tsuboi, 1959; Bullard and Coopery 1948; 
Tomoda and Ati , 1955; Dean , 1958] and polynomials [Oldham and Sutherland , 1955; 
Brown , 1956; Grant , 1957; Krumbein , 1959; Grant and Elsaharty, 1962; Mandel - 
baum 7 1963]. However, most of these have been devoted to the problem of ‘con- 
tinuation , [Peters } 1949; Hirvonen , 1952; Tengstrom, 1959; Orlin, 1959; Strakhov, 
1962] (extrapolation upward, to match airborne measurements, or downward, to 
determine crustal densities) and assume a fairly good distribution of observations. 
The most extensive application of Fourier methods has been by Kivioja [1962], 
who applied them to gravity anomalies over areas up to 10° by 35°. He found 
distance over which Fourier predictions could be extrapolated to be less than 4°. 

The most important application of gravity anomalies has been, and probably 
will continue to be, the determination of crustal structure in conjunction with 
seismic and geologic data [e.g., Woollard, 1959; Woollard et al, 1960; Press , 1960; 
Talwani et al, 1961; Oliver et al, 1961]. This circumstance suggests that the geology 
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and seismology of an area can be used to predict gravity anomalies where gravi- 
metry is lacking. These geologic methods have been applied most extensively by 
Durbin [1961] and Woollard [1962]. For the south central United States, where the 
rms anomaly is ±21 mgal, Durbin [1961] reports startlingly good results: the rms 
error of such predictions is only ±10 mgal. The situation does not appear to be 
so clear-cut in the more voluminous study by Woollard [1962], particularly for less 
stable regions. Also, one suspects that in many areas lacking gravimetry the 
geologic mapping w T ill also be inadequate. 

Statistical analysis has been applied to local variations of gravity anomalies 
by de Graaf-Hunter [1935], Hirvonen [1956, 1962], Kaula [1957, 1959a, 6], Baussus 
[1960, 1961], Moritz [1962a, 6, c], and Rapp [1962]. Fundamental to these analyses 
is the covariance function, as defined by (31). The principal defect of these analyses 
is the lack of adequate samples (as well as an adequate sampling theory) to ob- 
tain numerical estimates of covariance. The largest sample yet analyzed was by 
Kaula [1959a, b], consisting of nine areas each about 220 krn square in size and 
containing 52 to 140 members. Taking the data in blocks, rather than long lines, 
severely limits the number of independent estimates, because of the high correla- 
tion that will exist between two product pairs which are approximately parallel. 
The results obtained are shown in Figure 2, together with those obtained by Kaula 
[1957] from profiles across the United States and Hirvonen [1962] and Rapp [1962] 
from samples in Ohio and Finland, which unfortunately appear to be areas of 
smaller than average anomalies. Hirvonen finds from his sample that the formula 
for covariance over distance d , 

K(d) = K 0 /(l + c 2 d 2 ) (74) 


where K 0 and c are arbitrary constants, fits quite well, which is convenient for 
analytical development of other functions. However, a negative exponential form 
K 0 exp { — bd} or lower exponent on d in (74) appears to fit better to the larger 
sample of Kaula [1959a]. 

More commonly used than the covariance in earlier studies is the mean square 
anomaly for a square of side length s, related to the covariance K d by [Hirvonen, 
1962]: 


G m 


2 



WK d dr 


(75) 


where 


r = d/s 

W = (2tt - 8r + 2r 2 )r 0 < r < 1 (76) 

W = (2tt - 4 - 2r + 8 Vr 2 - 1 - 8 tan' 1 Vr 2 ^l)r r > 1 (77) 


If a single observation at the center is taken as representative of a square of side 
length 5, it wfill have a mean square error of representation: 

E 2 = Go 2 - G, 2 (78) 


The rigorous method of extrapolating and interpolating anomalies is to use the 
linear regression coefficients as specified by (46) and (47) [ Baussus , 1960; Moritz , 
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Fig. 2. Estimates of local and regional covariance of gravity anomalies. 
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19625]. Moritz also generalizes the linear regression to include the topography and 
anomalies corrected for linear correlation by equation 73. Kaula [1959a, 5] applied 
extensively to 1° by 1° square means the method of Markov estimation in terms 
of gravity and topography, which estimates for the (i + l)th member of a series 
of step-by-step extrapolations 


E{Ag »+i | hi+if Ag<, h;} 



xP\Ag i+1 = x, h i+1 | hi, Agr;} dx 
P{h i+1 | hi, Agf<} 


(79) 


As Baussus [1960] pointed out, this Markov estimation assumes that the co- 
variance function K(d) can be represented as a negative exponential of d . If K(d) 
is better represented by Hirvonen’s form (equation 74), appreciable improvement 
is obtainable when we use suitably spaced ‘next-to-nearest’ as well as nearest 
neighbors in estimation [Moritz, 19625]. 

The statistics of other functions can be derived as linear transformations of 
those for gravity anomalies, either with or without the intermediary of a spectral 
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representation: anomalies at higher altitude [Kaula, 1959a; Hirvonen , 1962; Moritz , 
1962a, c], geoid heights [KauZa, 1959a], and deflections of the vertical [Kaula, 
1959a; Kaula and Fischer , 1959]. 

Worldwide treatment . The statistical analysis of local variations described 
above has been extended to larger areas by Hirvonen [1956] and Kaula [1959a, b]. 
For covariance estimates, Kaula used, in addition to the local samples, eight regional 
samples covering an area about 10° by 10° with 56 to 115 members each (see 
Figure 2), plus a single world sample consisting of 569 5° by 5° mean anomalies 
based on observations in 18 per cent or more of their 1° by 1° squares. The co- 
variances obtained are shown in Figures 2 and 3. Table 1 gives the spherical 


TABLE 1. Degree Variances of Free Air Gravity Anomalies [Kaula, 19596] 
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harmonic degree variances a n 2 { A g } derived from the covariances by spectral analysis, 
as described by (32) and (29). These results are certainly of the correct order of 
magnitude, but the analysis needs to be redone with a properly randomized sample 
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[cf. Shapiro and Silverman , I960]. Hirvonen [1956] did not have the computer 
facilities necessary to make an extensive estimate of covariance, and so he used 
mainly mean square mean G 2 and error of representation E 2 , based on five regions 
averaging 8° by 8° plus the data of Tanni [1948]. His results imply a more gently 
varying field than Kaula’s. 

These authors and Cook [1950, 1951], Kaula [1957], and Molodenskiy et al. 
[1960] devote considerable attention to estimating the accuracy of determination 
of the geoid height and slope at specific points, given either the actual data distri- 
bution or a hypothetical dense net within a given distance of the point. Hirvonen’s 
and Kaula’s studies yield uncertainties that are probably too low, owing to the 
samples used and the assumption of randomness of 30° by 30° square means; 
Cook’s and Molodenskiy’s studies yield uncertainties that are probably too high 
because of use of the least-squares-determined spherical harmonic coefficients of 
Jeffreys [1943] and Zhonglovich [1952], respectively. 

The principal point of difference in various attempts to determine the external 
field from the incomplete information shown in Figure 1 has been on whether or 
not, and how, to utilize the topography as a means of interpolating and extra- 
polating over great distances. Dubovskiy [M olodenskiy , 1945], and Tanni [1948] 
used the topography by assuming zero isostatic anomaly in the unsurveyed areas. 
Kaula [19595], using the topography by applying the Markov extrapolation tech- 
nique of equation 79 to 5° by 5° means worldwide, probably obtained results 
similar to isostatic assumption although less susceptible to bias because of its 
broader statistical basis. In applying numerical integration, Heiskanen [1957] 
assumed zero free air anomaly with respect to the International Formula and 
Zhonglovitch [1952] assumed zero anomaly with respect to his least-squares-deter- 
mined third degree figure. Analyses limited to best fits to the observed gravimetry 
are those of Jeffreys [1943], Zhonglovich [1952] for harmonic degrees two, three, and 
four, and Uotila [1962]. All these solutions employed simple least squares — i.e., the 
harmonic coefficients were obtained as parameters z in (39) with coefficients C and 
covariance matrix W of the diagonal form specified by (44). To overcome the ill- 
conditioning due to neglect of covariance, Jeffreys [1943] grouped his data heavily 
into 30° by 30° square means and omitted all harmonic coefficients that the normal 
equation diagonal coefficient and constant indicated as making a small contribution 
before determining the remaining coefficients through degree 3 by least squares. 
Zhonglovich [1952] used the 10° by 10° square means and tried several least-squares 
solutions for different harmonics up through degree 4 before choosing as best a step- 
wise solution in which the coefficients for each successive degree 2, 3, 4 were deter- 
mined from the residuals with respect to previously determined degrees. Uotila 
[1962], benefiting from considerably more data, determined harmonic coefficients up 
through degree 4 by direct least-squares fit to 5° by 5° mean anomalies, holding 
the harmonics (n, m) = (1, 0), (1, 1), and (2, 1) zero and (2, 0), (3, 0), (4, 0) fixed 
at satellite-determined values. 

Jeffreys [1959, 1961] rejects use of the topography for long-range estimation; 
it is not clear whether he does so because the earlier determinations of gravity at 
sea obtained mostly positive anomalies, because the long-wave correlation is small 
(i.e., <r n (jh) small for small n in equation 64), because the applications of the topo- 
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graphy have involved the questionable isostatic or Markovian assumptions, or 
because, as a geophysicist, he is interested in estimating the amplitude, rather 
than the phase, of the variations. In any case it seems clear that (1) a solution 
optimum in the sense of quadratic sum minimization will be the one using the most 
information — i.e., the estimate by (72) will be better than that by (35); (2) the 
optimum estimates will be smaller in the mean square than the true values as a 
consequence of the smoothing effect of any prediction procedure; and (3) the ap- 
propriate measures of amplitude alone are the degree variances cr n 2 . 

It is difficult to characterize any of the existing solutions as approximations 
to the solutions described by (35) and (66), so that it is not clear to what extent 
they fall short of fully exploiting the available data. In principle, the general solu- 
tion up to about 12th degree harmonics should be easy with modem computing 
facilities: the half wavelength for the 12th degree is 15°, there are about 160 
15° by 15° squares with observations, and the 160 by 160 covariance matrix can 
be stored all at once in the core and inverted in a couple of minutes in an IBM 7094 
computer. The topography can also be incorporated in a more effective manner 
than before by using a new development thereof in spherical harmonics by Bruins 
[ V ening-M einesz, 1959] . 

USE OF SATELLITE ORBITS 

Observing system . The principal types of satellite tracking that have produced 
results useful for geodesy are essentially those described in the review by Kaula 
[1962a]: the 500-mm focal length f/1 Baker-Nunn tracking cameras; the 1000-mm 
focal length, // 5, modified aerial reconnaissance cameras, fixed and equatorially 
mounted; the Minitrack 108- and 136-Mc/s radio interferometers; and the Transit 
324- to 162-Mc/s radio Doppler trackers. The most important recent instrumental 
development has been the Anna geodetic satellite [Maconiber f 1963] which incorpo- 
rates magnetically oriented xenon-gas-discharge lamps capable of generating a 
flash of 8800 candle seconds along the axis, in addition to radio transponders and 
beacons. 

The distribution of the tracking stations of some of the principal systems used 
for geodesy is shown in Figure 4. The nonuniformity of distribution gives rise to 
some statistical problems because of the limited coverage of an orbit by a limited 
number of stations. Solving the spherical triangle formed by the orbital plane, 
the equator, and the meridian of a particular longitude X, we obtain 

os + / = sin -1 (sin <£/sin i) (80) 

Q — 6 = X — cot -1 (dz y/ tan 2 i esc 2 <j> — sec 2 i) (81) 

where i, to, /, and SI are, respectively, the inclination, argument of perigee, true 
anomaly, and longitude of the node of the orbital ellipse, and 0 is the Greenwich 
sidereal time. Since the inclination i varies but slightly, a station of position (<£, X) 
can observe the satellite only near certain values of co + / and ft — 0. For fairly 
well-distributed networks of about twelve stations, such as the Smithsonian 
Astrophysical Observatory Baker-Nu nn cameras and the Applied Physics Labo- 
ratory transit Doppler stations, the distortion due to this effect is probably slight 
except for series of camera observations dependent on solar illumination over dura- 



1LLIAM M. KAULA 



Some satellite tracking systems used for geodesy: C, Baker-Nunn cameras; D, transit Doppler receivers; I , Minitrack radio interferometers 

and Mots cameras. 
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tions short compared to a full cycle of nodal motion with respect to the sun — i.e., 
less than a few months. 

For the purpose of determining orbital variations, the random errors are quite 
satisfactorily small for all the systems mentioned above except the Minitrack 
interferometers, for which the effect of ionospheric refraction irregularities on ob- 
served directions are of the order of ±0.0005 to ±0.001, or ±100" to ±200". 
Ionospheric refraction is also probably a source of some systematic error in the 
Minitrack system. In the transit Doppler system, ionospheric refraction is believed 
to be largely eliminated by simultaneous transmission on two frequencies, and the 
most probable source of systematic error is shifting of the reference frequency pro- 
vided by the oscillator in the satellite, to the extent that the frequency is often 
taken as a separate unknown for each observed pass. The principal likely systematic 
error in camera observations is in the timing of observations of sunlit satellites, 
caused by difficulties of synchronization of the camera shutter with the clock and 
by error in the propagation delay time, possibly of the order of 0.001 sec. 

From the point of view of determining the gravitational field, the greatest 
systematic error common to all observational systems is error in tracking station 
positions. Error of positions of stations within the same triangulation system with 
respect to each other should be less than 20 or 30 meters; the principal possible 
exceptions are the positions of stations in South America and South Africa with 
respect to those in the same system in the northern hemisphere. Errors of position 
of triangulation systems with respect to each other and to the center of mass of 
the earth should be less than 100 meters for continental systems and less than 500 
meters for isolated stations. In addition to these errors that should exist, there 
have been cases where mistakes of as much as a kilometer have existed in station 
positions. 

Dynamics . The aspects of celestial mechanics important to understanding 
of close satellite orbits are explained in recent texts such as those of Baker and 
Makemson [1960] and Brouwer and Clemence [1961]. Satellite orbit dynamics with 
emphasis on geodetic applications is discussed by Kaula [1962] and Mueller [1963]. 

The theoretical problem of satellite orbits is to solve the equations of motion: 

f = V(F + fl.) + VA (82) 

where V is the gradient with respect to position, V, is the gradient with respect 
to velocity, F is the earth's gravitational potential as given by (9), R, is the gravi- 
tational plus radiation pressure potentials of third bodies (sun, moon, etc.), and 
R d is the atmospheric drag potential, a function of both position and velocity of 
the satellite. The information we wish to extract from satellite orbits rests in F, 
and so for the moment we neglect the other two potentials. The three second- 
order equations (82) are generally reduced to six first-order equations by change 
of variables: 


* = Z <?„(■) “7 (83) 

The simplest set of six variables s would, of course, be the position vector com- 
ponents { x , y, zj and the velocity vector components {x, y , z), referred to inertial 
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space with origin at the earths center. These six variables can be transformed to the 
six parameters of a Kepler ellipse with one focus at the origin: {a, e, i , £2, co, /}. 
The relationships between these parameters and the earth-fixed coordinates 
{ u , v, w) are shown in Figure 5. In the angle co + /, is the argument of perigee, 
the angle from the ascending node 12 to perigee, the point of closest approach of 
the ellipse to the origin, and / is the true anomaly, the angle from perigee to the 
satellite. Alternate ways of expressing the anomaly of the satellite are [Kaula y 
1962, p. 194] the eccentric anomaly E : 

tan \E = [(1 - «)/(l + e)] 1/2 tan \] (84) 

and the mean anomaly M : 

M — E — e sin E (85) 

The form of (83) in which s is the six Keplerian elements is advantageous 
because in the case of a purely central field, i.e., V is kM/r , the only nonzero rate 
of change is that of the anomaly; furthermore, for the mean anomaly M it is 
constant (Kepler’s third law) : 

n = - (kM) l/2 a 3/2 (86) 

where M on the right is the mass of the central body. 

In the case of a satellite moving in the actual potential field of the earth, by (9), 
for which the departures from a central field are 0(1(T 3 ), the Keplerian elements 
are still convenient because the rates of variation in (83) will be small except 
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for M , and hence representable as perturbations of the elliptic motion. The form 
of (83) for Keplerian elements is [Brouwer and Clemence, 1961, p. 289; Kaula, 
1962, p. 198]: 

a = (2/na)-(dV/dM) 

■ = 1 ~ e 2 dV _ (1 - e 2 ) 172 dV 
na 2 e dM na 2 e du 


cos i dV (1 - e 2 ) 1/2 dV 

na 2 (l — e 2 ) l/2 sin i di na 2 e de 


cos i 


dV 


dV 


di/dt na 2 (\ — e 2 ) l/2 s j n j na{\ — e 2 ) 1/2 sin i dil 


tl = 


1 


dV 


na 2 { 1 — e 2 )" 2 sin i di 

vr _ n H 2l. §¥. 

na e de na da 


(87) 


As is indicated by Figure 5, it is a purely geometrical problem to transform 
V (r, </>, X) as given by (9) to a form V (fl, i, oj, a, e , M ) suitable for taking deriva- 
tives in (83). The solution which is given by Kaula [1961a]: 

V = kM \l + l E (f*)’ E N nm E F nmv {i) • E G n „(e)S nm „( u>, M, 0, 0)1 (88) 

a „ — 2 / n» " 0 p-0 « a, J 

where the normalization factor 


iV n . 


_ (2ra + l)(n - m)l (2 - 3 0m ) 1 
(n + m)! 


f w (o = E 


(2n - 2Q! 


T - t\ (n - t )\ (n - m - 202 2 "' 


_ * »— m—21 • 

■ sin l 


cos 


‘i Z 


n — m — 2t + s 
c 


m — s 
Ip — t — c 


(-i )• 


(89) 


(90) 


in which t is summed from 0 to p or k (defined after equation 11), whichever is 
less; s, from 0 to m; and c , over all values making the two binomial coefficients 
nonzero; 


G np( ^ n) (e) = (1 - 6 2 ) 


p'-i 

2^ (1/2) — n 

rf-0 


71—1 
2d + n — 2p'J 


2d + n - 2p' 

d 


2d+n— 2p' 


(91) 


where p' = p, p < n/2, p’ — n — p, p < n/2, for q = 2p — n;ior q ^ 2p — n, 
more complicated infinite series are required [Kaula, 1961a, equations 24-26]; and 


S„ mvq ( to, M, 2, 0) 

c nm 

— Snm 

s., m 


cos {(n — 2)co + (n — 2p + q)M + m(Sl — 0)} 


+ 


(»— m ) odd 
(n— m) even 


L C n 


sin {(n — 2p)a> + (n — 2p + q)M + m(Q — 0)} (92) 


(«— m) odd 
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Other expressions are given by Musen [1960] and Groves [I960]. If V as expressed 
by equations 88-92 is differentiated and the derivatives are placed in (87), the 
form of S nmvQ in (92) indicates that a, e, di/dt must be sinusoidal. If we assume 
that these sinusoidal variations are small enough that a, e, i can be considered 
as constant on the right of (87), then (92) further indicates that only «, Q, and M 
have constant terms, and that these terms arise from cases where n is even, m 
is zero, p is n/2, and q is zero. Hence to a linear approximation, (87) can be inte- 
grated assuming a, e } i, n constant (except for n in M , the variation of which is 
obtainable from equation 86 and a) and M, w, 12 secularly changing. 

However, since C 2t0 is 0(10“ 3 ) while all the other C nm and S nm are 0(10“ 6 ), 
a linear approximation does not suffice for C 2t0 : terms with coefficient C 2 , o are 
required. As usual when a problem becomes nonlinear, it becomes appreciably 
more complicated, and care must be taken in defining the constants of integration. 
The obvious choices are the values of the Keplerian elements at a particular instant 
of time, i.e. the osculating elements. The osculating elements are the constants 
of integration when the problem is solved in a purely numerical manner. 

When the problem is developed further for an iterative or analytic solution, 
however, it is usually found more convenient to define the constants as elements 
of a fictitious reference orbit or an intermediate orbit. In some theories this inter- 
mediary is defined geometrically, as, e.g., representing all constant and secular 
terms of the actual orbit; in other theories it can be defined dynamically as cor- 
responding to a constant part of the potential V. Other aspects in which theories 
may differ are the coordinate system; the independent variable; the orbital ele- 
ments— often a canonical set is used, i.e. one such that the C if - in (83) consists of 
only one nonzero constant per row; and the point at which the development of 
the problem changes from algebraic to numerical. This large variety of possibilities 
gave rise to a large number of theoretical papers in the period 1957-1961 on the 
orbit of a close satellite of an oblate planet to at least 0(J 2 2 ) in secular motions 
(J 2 is — V5 C 2 , o). The theories that probably have been used the most are those 
of Musen [1959], who adapted the Hansen lunar theory to a form suitable for solu- 
tion by iteration; Brouwer [1959], who applied Von Zeipeks method of canonical 
transformation with a purely Keplerian intermediary; Kozai [1959a, 5], who ex- 
tended the Lagrangian equations 87 to higher-order terms; V inti [1959, 1961] 
and Izsak [1960], who both separated the equations of motion by using ellipsoidal 
coordinates; King-Hele [1958], who employed a Keplerian ellipse of fixed inclination 
and perigee argument as intermediary and solved in successive approximations 
according to powers of J 2 and e; and Merson [1961], who made a development 
similar to King-Hele's starting from osculating elements when the satellite is 
at the node. (See Kaula [1962] for more description and comparison.) 

Since 1960, more theoretical attention has turned to the problems of resonance 
associated with close satellites. The problem of critical inclination in the vicinity 
of sec" 1 •%/ 5, or 63°26', at which perigee motion is zero, is of relatively little geodetic 
interest; the most complete solution is probably that of Izsak [1962], as extended 
by Aoki [1963]. Of much more geodetic significance are the 24-hour orbits (semi- 
major axis 42,000 km), which resonate with tesseral harmonics C nmj S nm for which 
n — m is even, since 


f2 + ci> + M— 0£^O 


(93) 
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The problem of orbits to which (93) applies with potential 

y = ~~ + R = ^ 22 F 220 G' 200 S 220 o] (94) 

has been analyzed by Blitzer et at. [1962, 1963], who use a linearized system of 
differential equations; by Musen and Bailie [1962], who isolate (Q + a> + M — 6) 
as a canonical element and dev elop the Hamiltonian in powers of the ratio 
(< dR/dL)/(d 2 R/dL 2 ), where L = \/kMa; and by Morando [1963], who applies Von 
ZeipePs method and develops the determining function in powers of (C 22 2 + & 22 2 ) 1/4 • 

Elaborately developed theories such as those of Musen [1959] and Brouwer 
[1959] are advantageous to analyze secular changes for zonal harmonics and to 
conserve computer time when observations are infrequent (say, an average of 
less than one per hour). Computation of orbital arcs of a few days or less using 
frequent radio tracking is still done mostly by numerical integration. 

As necessary as a correct dynamical theory, and often more laborious, are 
differential correction systems to determine from observations the numerical 
values to use in the theories. Examples of differential correction programs for close 
satellites are those of Veis and Moore [1960] and Merson [1963]. The geometrical 
aspects of differential correction schemes are fairly complicated but straightforward 
(see Veis [1960] and Kaula [1961a, 1962] for details). The statistical aspects are 
not so complicated, but are much less satisfactorily treated, partly because a 
rigorous treatment would require excessive computer storage and time, and partly 
because adequate statistical models are lacking. Satellite orbit observations always 
have an appreciable amount of serial correlation, principally because of the ina- 
bility to account for air drag. However, the statistical analysis of satellite orbit 
accelerations has received relatively little attention, the only active worker thereon 
being Moe [1963]. No differential correction method takes into account correlation 
between observations of different passes, as suggested in the section of this review 
on time series. The most elaborate treatment existing [Kaula, 19635] allows for 
correlation within a pass and applies various other devices such as giving greater 
weight to the across-track than to the along-track component of an observation; 
preassigning variances and covariances of parameters being determined in the 
analysis, so that mathematically they become 'observations’; and weighting ob- 
servations according to their distribution with respect to phase angles believed 
critical. 

Determination of zonal harmonics. For the principal secular or long-period 
effect of a zonal harmonic, m = 0, the disturbing terms in (88) can be taken as: 

Long period R„ = -kMJ n ^ •F„. t (t)G„ t(2i _.,(e)[ 1 ’ _ n 6Ven (95) 

( 2 sin o), n odd 

where J n is — a/2 n + 1 C n , 0 . The linear perturbations can be obtained by using 
this disturbing function in the Lagrangian equations of motion 87; in addition, 
nonlinear terms J 2 2 and «/ 2 J„, n odd, need to be taken into account. 

The customary method of determining the JJs from satellite orbits is to 
analyze the long-term variation of orbital elements determined by a differential 
correction fit to observations over a few days at a time. Precautions that must be 
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taken in analyzing secular changes of the node 12 and perigee co to determine even 
degree zonal harmonics include: 

1 . The set of satellite orbits used should have a variety of inclinations suffi- 
cient to separate the different harmonics. 

2 . The constants of integration determined by analyzing observations must 
be consistent with the algebraic form of the terms containing J 2 ( C 2( o 2 in our 
notation). 

3. The mean value of the constants of integration must be accurately de- 
termined for the duration from which the secular rates 12 , co are determined. It is 
also to be emphasized that mean values of the constants of integration will differ 
from mean values of osculating elements for some theories (e.g., the inclination 
in the theory of King-Hele [1958], or the mean motion in the theory of Brouwer 
[1959]). Averaging of the elements a and e is important in order to remove secular 
drag effects; if the a, e used are not average values, but rather values for some 
epoch noncentral to the duration used, about the best that can be done is to use 
the perturbation in mean anomaly to correct node and perigee assuming the 
perigee height to remain fixed [O'Keefe et al, 1959; Kozai , 1962, 1963a]: 

< 96) 

The effects of errors Ac, A i, A a in the orbital elements on determination of 
the secular rate of the node are [Kozai, 1962]: 

. T 4e 7 1 

A12 ~ ^ y ~ — 2 Ac — tan i A i — — Aa (97) 

For some satellites of high inclination, error in determining zonal harmonics will 
thus come more from error in the mean inclination than in 0 itself. 

4. If luni-solar attraction, radiation pressure, and other perturbations are 
not removed in determining the mean values of the constants of integration, they 
can distort determination of the rates 12 and co not only through purely secular 
effects but also through periodic perturbations. A periodic perturbation A ( 12 , co) sin 
[ Kt — X ) will affect the apparent secular rate from observations lasting from U 
to t 2 by an amount 6 ( 12 , A): 

5 ( 6 , *) = M ■ - sjn ,_ U (98) 

t 2 t\ 

5. If the perturbations are removed in determining the constants of inte- 
gration, in addition to direct effects Aj(12, co), the interaction of perturbations Ae, 
At with the secular effect of J 2 may cause an indirect effect A 2 (12, co) large enough 
that it should be taken into account: 

A 2 (fl, u) = f A edt + ^ f ~ f Ai dt (99) 

6 . Satellites with low perigee, nonspherical shape, and large area-to-mass 
ratio should be avoided because of the difficulty in calculating drag and radiation 
pressure effects on the mean elements and the secular rates, particularly on co. 
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7. Use of orbits provided by routine prediction services should be avoided, 
because of possible imprecision in the definition of orbital elements and inaccuracy 
in the determination thereof. 

The most important recent analyses of secular motions have been by Kozai 
[1962] and King-Hele et al. [1963]. Kozai uses Cl and ci of thirteen satellites ranging 
from 32.9° to 66.8° in inclination. However, he weights the data inversely pro- 
portional to the squares of the standard deviations, which vary greatly, so that 
the result is almost entirely determined by the Cl of 1960i 2 (inclination 47.2°) and 
the Cl and d> of 1959c*! and 195977 (inclinations 32.9°, 33.4°). King-Hele et al. [1963] 
used only the nodal motion Cl of seven satellites at widely spaced intervals of in- 
clination from 32.9° to 97.4° and weighted each satellite about equally. However, 
to obtain inclinations above 53.8°, they used elements provided by routine pre- 
diction services. The disagreement between the results of Kozai [1962] and King- 
Hele et al. [1963] in Table 2 at this late date is disappointing. What are required 
are accurately determined orbits of high inclination. A limited amount of the 
requisite data are now available in the form of 552 precisely reduced Baker-Nunn 
camera observations over 54 days of 1961aS a , which has inclination 95.9°, perigee 
height 3500 km. The observed nodal motion Cl determined as a byproduct of 
analysis for tesseral harmonics [. Kaula , 1963d] is +0.21037° db 0.00010/day. 
Including a luni-solar effect of — 0.00006 °/day, the calculated motion Cl is 
+0.21034 °/day using the J 2 through J s of Kozai [1962], and 0.21056°/day using 
the J 2 through J 12 of King-Hele et al. [1963]. 

The determination of the odd degree zonal harmonics is somewhat easier, 
since there do not appear to be any other significant effects of period 2 tt / u >. The 
principal precautions are to include the J 2 J n terms, following (99), and the drag 
correction of (96) [Kozai, 1959a]. Results by the principal investigators are given 
in Table 2. 


TABLE 2. Zonal Harmonic Coefficients of the Gravitational Field 


Coefficient* 

Newton et al. 
[1961] 

Smith 

[1961, 1963] 

Kozai 

[1962] 

Shelkey 

[1962] 

King-Hele et al. 
[1963] 

Anderle and 
Oesterwinter 
[1963] 

Ji X 10* 


1083.15 

1082.48 

1082.61 

1082.86 

1082.47 

J, X 10* 

—2.36 

-2.44 

-2.56 

-1.94 


-2.48 

Ji X 10* 


-1.4 

-1.84 

-1.52 

-1.03 

-1.40 

Ji X 10* 

-0.19 

-0.18 

-0.06 

-0.41 


-0.14 

Ji X 10* 


0.7 

0.39 

0.73 

0.72 


Ji X 10* 

-0.28 

0.30 

-0.47 




Jt X 10* 



-0.02 


0.34 


J, X 10* 



0.12 




Ju X 10* 





0.50 


Ju X 10* 





0.44 



*J n - “V2 n + 1 C n# 


Determination of GM. The rapid motion of artificial satellites suggests that 
their mean motion may serve as an accurate method of determining GM, the product 
of the gravitational constant and the earth's mass, through Kepler's equation 86. 
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Drag and radiation pressure have their greatest effects on the mean anomaly, 
however, and so it is particularly desirable to have a satellite for which these effects 
are minimized and readily calculable: spherical, low area-to-mass ratio, small 
eccentricity, high inclination, perigee height above 1000 km. Determinations of 
GM from such satellites have been made as byproducts of tesseral harmonic 
analyses by Kaula [19635, d], using camera observations, and by Anderle and 
Oesterwinter [1964], using Doppler observations. The determination from camera 
data depends on the geodetic triangulation connecting tracking stations to estab- 
lish scale; that from Doppler data may be affected by ionospheric refraction effects. 
A more distant satellite would reduce these effects: in this category are the use of a 
lunar probe [Hamilton et ah, 1963] and those methods measuring the distance of 
the moon by radar [Yaplee et ah, 1963] and by triangulation [Fischer, 1962]. The 
radar method now obtains an internal accuracy of ±200 meters for the distance 
of the moon by correcting for the variations in topography on the moon’s surface. 
Both the radar and triangulation methods are affected by error in the radius in the 
moon in the direction of the earth, for which the scatter of different determinations 
is about 2 km, equivalent to a variation of 0.00007 X 10 14 m 3 /sec 2 in GM. Another 
possible source of error is the ratio p of the moon’s mass to the earth’s, appearing 
in the modified Kepler equation 

GM = nV i 1 ; (100) 

(1 + m) 

where 0 is the effect of the sun on the mean distance. When the astronomical unit 
found by radar measurements to Venus is used, the values of /T 1 deduced from the 


TABLE 3. Determinations of GM 


Method 

Reference 

Sources of Error 

GM 

10 20 cm 3 sec -2 

Terrestrial geodesy 

Fischer [1962] 

fTrianguIationl 

3.986040 


Kaula [19616] 

[Gravimetry f 

3.986020 ± 0.000028 


a e from Kaula [19616] 


3.986043 


y e from Uotila [1962] 



Lunar motion and 

Yaplee et at. [1963]* 

Lunar radius 

3.986057 

radar distance 




Lunar motion and 

Fischer [1962]* 

/Lunar radius \ 


triangulated 

Crommelin a 

(Triangulation J 

3.986451 

distance 

Fischer [1962]* 


3.986078 


O’Keefe and Anderson a 



Lunar probe and 

Hamilton et al. [1963] 

Observational 

3.986016 =fc 0.000025 

Doppler 


Station position 


Close satellite and 

Anderle and Oesterwinter 

Observational 

3.985889 

Doppler 

[1963] 

Orbit perturbations 


Close satellite and 

Kaula [19636], 1960 1 , 

JTriangulation \ 

3.986037 ± 0.000012 

camera 

Kaula [1963d], 196 la % 

(Orbit perturbations/ 

3.985993 ± 0.000011 


* These references are sources for only the distance to the nearest point of the moon; a 
lunar radius of 1738.7 km and a lunar mass of 1/81.3015 have been used to calculate the values 
of GM given in the last column. 
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lunar inequality (the monthly revolution of the earth about the center of mass of 
the earth and moon) based on observations of Eros range from 81.26 to 81.36. 
A new determination of the inequality from Doppler observations of the Mariner 2 
Venus probe yields a preliminary value of 81.3015 ± 0.0033 [Hamilton et aL , 1963]. 

Determination of tesseral harmonics. The determination of tesseral harmonics 
from satellite orbits depends principally on the terms in the disturbing function, 
(88), for which the argument is {(n — 2p)co + m(S2 — 6)\ y and for which \q\ is 
small. Because of the unavoidable inclusion of the earth's rotation rate 6, the 
frequencies involved differ by orders of magnitude from those involved in analysis 
for zonal harmonics, and so the problems are quite different. On the one hand, 
long-term-drag, luni-solar, and radiation pressure effects are of little influence; 
on the other hand, errors in station position and nonuniform distribution of ob- 
servations become important. Usually the frequency of observation is not much 
higher than the frequencies of the orbital perturbations caused by the tesseral 
harmonics; hence it is not possible to employ the method of first determining 
smoothed osculating elements from the observations and then analyzing variations 
in these elements. The analysis must be made directly in terms of the observations 
themselves; i.e., there must be formed a partial derivative of each observation with 
respect to each of the parameters sought. A rigorous solution will then follow that 
of the times series of equations 57-60, in which x f (t) are the Keplerian elements, 
p t are the gravitational coefficients, ?/ 7 - are the observations, and q m are the cor" 
rections to station position coordinates. The difficulty in the rigorous solution, as 
stated, is taking into account the covariance between observations, since this 
entails arrays of dimension comparable to the number of observations. Further- 
more, the statistical analysis of drag effects has not been developed sufficiently 
[Kaula, 1961c; Moe , 1963]. Consequently, in practice all analyses have neglected 
covariance between observations of different passes and have applied various de- 
vices to minimize the effect of this neglect: 

1. Higher weighting of the across-track than of the along-track component 
of an observation is used, because a proportionately much greater part of the drag 
effect is in the mean anomaly than is true for the gravitational effects [. Kaula , 
19635, Izsak , 1962]. 

2. The duration for which a set of reference elements, or constants of inte- 
gration, are determined is limited to one to four weeks. This measure is also de- 
sirable to keep the residuals down to not more than a small multiple of the antici- 
pated gravitational effects. 

3. Arbitrary polynomials are used to represent some of the variation in mean 
anomaly. 

4. The observations are weighted inversely to their density with respect 
to the angle (O — 8) in an attempt to restore some of the orthogonality of gravi- 
tationally caused variations to the variations caused by drag [Kaula, 19635]. 

The nonuniform distribution of observations arises from the geometrical effect 
described by (75) plus, in the case of camera observations, from dependence on 
clear weather and solar illumination of the satellites. This nonuniform distribution, 
the large number of parameters involved, and the similarity of the perturbations 



538 


WILLIAM M. KAULA 


by different gravitational terms (l, m) and (n, m) for which (l — n) is even, all 
combine to produce an ill-conditioned least-squares solution from observations 
of a single satellite. The desirable remedy is to combine several arcs of several 
satellites of widely varying inclination in one solution. However, such a solution 
would involve well over 100 unknowns, so further compromises must be made. 
These compromises have been of two principal types: 

1. Stepwise least squares. This method is employed by Izsak [1963], Anderle 
and Oesterwinter [1964], and others. For camera observations there are two steps: 
in the first step, the reference elements for each orbital arc are determined holding 
fixed the gravitational coefficients and station positions; in the second step, the 
orbital elements are held fixed, while the gravitational coefficients and station 
positions are determined from the residuals of the observations with respect to 
these reference orbits. For Doppler observations, there are three steps: the above 
two steps are preceded by a step in which frequency constant and frequency drift 
corrections are determined for each pass. The likely defect of these methods is 
that the parameters determined in the earlier steps will absorb some of the effects 
of the gravitational coefficients and station position shifts. This defect may be 
particularly severe in the case of orbital arcs with relatively few camera observa- 
tions. Izsak [1963] omits the worst of such arcs, but this measure introduces further 
chance of bias because the arcs so omitted will tend to be those for which the apogee 
is in the southern hemisphere, away from most of the stations. Both Izsak [1963] 
and Anderle and Oesterwinter [1964] make solutions in which all station positions 
are allowed to move freely with respect to each other; however, computer program 
limitations do not permit them to determine all the gravitational coefficients 
compatible with the accuracy implied by this method. 

2. Preassigned covariance matrix [ Kaula , 19636, d]. For each orbital arc, the 
gravitational coefficients and station positions are determined at the same time 
as the reference elements. To keep the solution from ‘blowing up’ because of ill- 
conditioning, a covariance matrix W, is preassigned to the gravity and position 
parameters in accordance with (50). The variances in W, for the gravitational 
coefficients are based on the degree variances in Table 1; those for the positions of 
major datums are based on the results of the world geodetic system solution of 
Kaula [19616], and those for the positions of isolated stations on statistical analyses 
of deflections of the vertical. Stations connected by geodetic triangulation are 
assumed to translate together, consistent with an accuracy of about d=20 meters. 
The principal defect of the preassigned variance and covariance method is that it 
may influence the relative magnitude of the results, particularly for harmonics 
causing similar frequency perturbations in the orbit. 

Earlier estimates of tesseral harmonics from satellite orbits [Izsak, 1961a; 
Kaula , 1961c; Kozai , 1961; Newton, 1962] yielded a wide scatter of results and are 
of interest now only for some discussion of methods. Table 4 gives the most recent 
results of most of the principal investigators; those by Kozai [19636], Izsak [1963], 
and Kaula [1963d] are based on camera data; those by Anderle and Oesterwinter 
[1964] and Guier [1963] are based on Doppler data. In addition to the gravitational 
coefficients given in the table, Kaula [1963d] determined 10 tesseral coefficients 
of degrees 5 and 6, plus 18 datum coordinate shifts; Anderle and Oesterwinter 
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[1964] also determined 54 station coordinate shifts. The most recent determinations 
from terrestrial data [Kaula, 19615; Uotila , 1962] are also given in Table 4. Figure 6 
is the geoid map corresponding to the solution by Kaula [1963d]. 


TABLE 4. Harmonic Coefficients of the Gravitational Field 


Coefficient 

Kaula 

[1961b] 

Uotila 

[1962] 

Kozai 

[19636] 

Izsak 

[1963] 

Anderle and 
Oesterwinter 
[1964] 

Kaula 

1963d] 

Guier 

[1963] 

C 20 X 10' 

-484.23 

(-484.10) 

(-484.10) 

(-484.10) 

-484 09 

-484.08 


C 22 X 10' 

0.75 

0.69 

1.11 

1.50 

2.85 

1.88 

2.60 

Sa X 10' 

-0.38 

-2.25 

-1.47 

-0.62 

-1 53 

-1.38 

-0.99 

C 30 X 10' 

0.78 

(0.91) 

(0.97) 

(0.97) 

0.94 

0.97 


Cn X 10* 

1.03 

0.10 

1.75 

1.04 


1.52 

1.64 

S,i x 10 ' 

0.39 

-0.63 

0.14 

0.06 


0.14 

0.18 

C H X 10' 

0.97 

1.19 

0.34 

0.27 


-0.02 

0.84 

S a X 10* 

-0.03 

-0.22 

-0.24 

-0.53 


0.42 

-0.07 

C 33 X 10' 

0.57 

1.50 

-0.45 

0.51 


0.70 

1.06 

S„ X 10' 

1.40 

2.44 

0.60 

0.89 


0.76 

1.01 

C.o X 10« 

0.25 

(0.57) 

(0.61) 

(0.61) 

0.47 

0.67 


C« X 10' 

-0.63 

-0.15 

-0.30 

-0.30 

-0.71 

-0.33 

-0.60 

X 10' 

-0.15 

-0.20 

-0.46 

-0.34 

-0.40 

0.37 

-0.49 

C« X 10' 

0.46 

0.82 

-0.18 

0.16 

0.45 

0.01 

0.27 

S« X 10' 

0.42 

0.46 

0.21 

0.55 

1.20 

0.35 

1.19 

C„ X 10' 

0.51 

1.20 

0.59 

0.36 

2.64 

0.17 

1.33 

S a X 10' 

-0.01 

-0.63 

0.02 

0.25 

-0.60 

0.41 

-0.05 

C u X 10« 

-0.10 

0.66 

0.78 

0.46 


-0.01 

-0.37 

Su X 10' 

0.36 

-0.12 

1.26 

0.77 


0.18 

0.31 


USE OF ASTROGEODESY 

The oldest sources of information about the variations in the gravitational 
field are the slopes of the geoid with respect to the surface of a reference ellipsoid 
determined from the differences between the astronomic and geodetic positions: 




( 101 ) 


r] = (\ a — X 0 ) COS <f> 

Along a line of geodetic triangulation or traverse in azimuth A the slope will be: 

t = r\ sin A + £ cos A (102) 

and for the change in geoid height over an arc distance 'F there is conventionally 
applied the simple integration 


N b = N, 


-*/' 

Jo 


tdl 


(103) 


where R is the radius of the earth. 

However, if the triangulation from which the geodetic positions (<t> 0f \ 0 ) 
were calculated with a scale correction R/(R + h) for reduction to sea level 
altitude h (the 'development method’), then an error will accumulate because 
of neglect of the further scale correction R/(R + N) for reduction to ellipsoid 
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Geoid heights referred to an ellipsoid of flattening 1/298.24, based on camera observations of satellites 1959ai, 1959ij, 1960i«, 19615j 

1961«Si [ Kaula , 1963d], 



THE EARTH’S GRAVITATIONAL FIELD 


541 


(the ‘projection method'). After a few thousand kilometers, this error will make 
a difference of several tens of meters in N B calculated by (103). The elimination 
of this discrepancy is known as the Molodenskiy correction; given t, N calculated 
by (102) and (103) from development-computed </>„, A fl , the correct slope 0 and 
height f will be [Molodenskiy et al., 1960, p. 33] 

1 C* 

0 B = t B + (fix — t A ) cos ^ + (f A — N A ) sin \j/ + ^ j N cos (l 
ie = N b + R(t A — e A ) sin yf/ + (f A — Nx) cos $ + [ N sin ( l 

Jo 

If geodetic control computed by the development method is adjusted forcing 
conditions on large circuits properly applicable only to projection-computed con- 
trol, then (104) will not remove all error. However, circuits large enough for this 
distortion to be significant are rare. 

The most extensive application of (101) through (104) to astrogeodetic data 
has been by Fischer [1959a, b; 1960a, b) 1961]. Figure 7 shows the distribution of the 
data she has applied. The limited geographic extent of the astrogeodetic data 
makes it primarily of value as an independent check on the results obtained by 
gravimetric and satellite means. The significance of the various spacings of the 
astronomic stations in Figure 7 is that the principal source of error in determining 
astrogeodetic geoid height differences is the error in determining the continuously 
varying geoid slope by interpolation between the observed point values. For sta- 
tion spacings S in kilometers, the mean square expected error in the geoid height 
E„ { e 2 } derived from the degree variances in Table 1 (extended to higher degrees) 
is approximated by [. Kaula , 19616]: 

E'{e x 2 ] « (0.0195 - 1.4) 2 (105) 

This rule probably gives an underestimate for the geoid height error resulting from 
arcs carried through rugged topography such as that down the west coast of South 
America. However, it seems safe to conclude that the rms error for the relative 
location with respect to each other of points in the major systems shown in Figure 7 
is of the order of ±15 meters in the radial coordinate, as well as the two hori- 
zontal. This estimate is confirmed by misclosures of less than 25 meters for 10,000- 
km loops around the Caribbean [ Fischer , 1959] and around the Black and Caspian 
seas and through Turkestan [Fischer, 1961]. (See Bomford [1960; 1962, pp. 143-159, 
325-327] for more detailed discussion of triangulation and geoid height accuracy, 
and Rice [1962] for an example of an optimum geoidal section survey.) 

There are several ways of combining astrogeodetic and gravimetric geoid 
data. In the USSR, gravimetry is used to interpolate deflections between astro- 
geodetic stations [Molodenskiy et al , 1960, pp. 125-141]. The traditionally advo- 
cated method of combination is to compare the astrogeodetic and gravimetric geoid 
heights at a few carefully selected points and then to minimize the discrepancy 
between the two by adjusting the datum position and ellipsoid parameters 
[Heiskanen and Vening-Meinesz , 1958, pp. 299-310]. This method has been applied 
extensively by Rice [1952] to sixteen stations in the United States and by 
Szabo [1962] to six stations in the United States and twenty-three stations in 


- i p)dl 

(104) 

- f)dl 
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Eurasia. However, since there is still appreciable error in the gravimetrically 
computed deflection at even the best points (probably at least ±1.5") the logical 
conclusion is to compare astrogeodetic and gravimetric geoids wherever the former 
is available. Zhonglovich [1956] made such a comparison by minimizing the sum 
{(la — Vo) 2 + (L — f 0 ) 2 } for the mean values of ninety-six 4° by 4° squares 
covering North America, using his geoid l Zhonglovich , 1952] for The mean 

curvature of the ellipsoid he believed best fitting yields an equatorial radius of 
6,378,104 ± 42 meters when the flattening of 1/298.3 is enforced. Fischer [1960, 
1961] minimized the sum (Na — N g ) 2 at 301 points at 5° intervals throughout 
the area pictured in Figure 7, using N g calculated by Kaula [1959a, b\. Fixing the 
flattening at 1/298.3, the equatorial radius obtained varied from 6,378,160 to 
6,378,166 meters on different assumptions. Kaula [19616] used the same data as 
Fischer, but applied a much more detailed statistical treatment, in which 106 
condition equations were written requiring that the corrected differences in astro- 
geodetie geoid height between mean values for 10° by 10° squares agree with 
the difference calculated from the corrected harmonic coefficients of the gravi- 
tational field. The quadratic sum minimized was 

minimum = x A W A ~ 1 x A + x G W 0 _1 x G + x.W. -1 *. (106) 

where x A comprises the corrections to the 106 astrogeodetic geoid height dif- 
ferences; W A is the associated covariance matrix, taking into account interpolation 
error, error of representation, and the discrepancy between the actual geoid and 
that represented by harmonics to the 8th degree; x G comprises the 81 spherical 
harmonic coefficients through the 8th degree; W<? is the covariance matrix produced 
by the analysis of gravimetry by Kaula [1959a, 6]; and x, and W, pertain to supple- 
mentary measurements of secular and long periodic satellite motions and geoid 
height matching between adjacent but unconnected datums. A measure of the 
agreement between astrogeodetic and gravimetric geoids was the value obtained 
for the quadratic sum of (106); it was 44 per cent higher than the mean x square 
expectancy. Increasing the standard deviations accordingly, 6,378,163 ± 15 meters 
and 1/298.24 ± 0.01 were obtained for the ellipsoid parameters. The values of the 
harmonic coefficients obtained through the 4th degree are given in Table 4. The 
standard deviations obtained for the tesseral harmonic coefficients of degree n 
averaged about ±0.9 X l0~ 6 /(n — !)• The rms discrepancies from the satellite 
solution of Kaula [1963d] are ±1.0 X 10“ 6 for degree 2, ±0.6 X 10 -6 for degree 3, 
and ±0.3 X 10'° for degree 4, so that the disagreement is one which can be reason- 
ably expected. 


CONCLUSIONS 

The application of more elaborate statistical techniques made possible by 
modem computers could very probably extract more information about the 
gravity field from existing data, both gravimetric and satellite, and would facilitate 
the planning of programs of additional observations. However, these elaborations 
would be added to what is already a rather complicated task of measurement, 
data processing, and theoretical analysis, such that both the benefits and penalties 
of any modifications are difficult to predict. Also, this review has not touched on 
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the geophysical application of the gravity field. It can be argued that an appreciable 
effort to determine the gravity field is not worth while because knowledge is lacking 
of the long-term rheology of the earth’s interior as well as the mathematical tech- 
niques to apply all but the simplest rheologies. Certainly much of what has been 
discussed in this review is not necessary for the application of gravimetry to study 
local and regional variations in the crust. What part of it is useful to study vari- 
ations on a continental or oceanic scale and the implications thereof as to the 
state of the mantle is perhaps worthy of consideration: To what harmonic degree 
should the field be developed to apply to ideas of a weak upper mantle or to con- 
vection currents? Are the degree variances a useful tool? Should coefficients be 
derived by techniques which yield mean square values smaller than those of the 
actual coefficients? The answer is probably the customary one that we cannot 
predict what the future will want, and hence should strive for the most accurate 
representation possible. 
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